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RESEARCH  RESULTS 


During  the  period  of  sponsorship  by  AFOSR  Grant  Number  AFOSR  F 49620-93- 1-0061  we 
have  undertaken  and  completed  a  variety  of  research  projects.  Here,  we  briefly  describe 
our  research  efforts.  Below,  we  will  give  details  concerning  four  of  these  projects. 

We  have  undertaken  a  systematic  analytical  and  computational  study  of  least- squares 
finite  element  methods  for  incompressible  viscous  flows.  Details  of  this  research  axe  given 
below  and  in  papers  (l]-(5j  listed  in  the  bibliobraphy.  We  have  developed  and  analyzed  the 
first  optimally  accurate  least-squares  finite  element  method  for  the  Navier-Stokes  equations 
with  velocity  boundary  conditions.  Our  analyses  include  the  nonlinear  case  ([2]  and  [5]) 
as  well  as  the  linear  Stokes  case  ([1],  (3j,  and  [4]).  In  this  respect,  we  have  performed 
perhaps  the  first  rigorous  analysis  of  a  least  squares  finite  element  method  for  a  nonlinear 
problem.  We  have  also  performed  extensive  computations  which  serve  to  demonstrate 
the  implementation  of  the  techniques,  as  well  as  to  provide  illustrations  of  the  theoretical 
results. 

A  second  major  undertaking  is  the  optimal  control  of  incompressible  viscous  flows.  De¬ 
tails  are  given  in  papers  [6],  [8],  [10],  [14],  [15],  and  [17]  listed  in  the  bibliography.  We  have 
considered  both  boundary  and  shape  controls.  We  have  developed  and  implemented  algo¬ 
rithms  and  analyzed  a  multidisciplinary  optimization  problem  involving  coupled  solid/fluid 
heat  transfer  [17].  We  have  an  on-going  project  in  the  development  and  implementation  of 
algorithms  for  the  shape  control  of  flows  [6].  We  have  extended  our  results  for  the  Navier- 
Stokes  case  to  a  model  for  incompressible  viscous  flows  due  to  Ladyzhenskaya  which  fea¬ 
tures  a  nonlinear  constitutive  law  [8].  One  highlight  of  our  research  is  that  we  have  given 
the  first  rigorous  analyses  for  a  shape  control  problem  for  the  nonlinear  Navier-Stokes 
problem;  we  provide  some  details  below  and  in  [14]  and  [15]. 


As  an  extension  of  our  work  on  optimal  control  of  incompressible  flows,  we  have  been 
able  to  develop  and  abstract  theory  for  the  finite  dimensional  approximation  of  a  class  of 
nonlinear  optimal  control  problems.  Details  are  given  below  and  in  the  paper  [13]  listed 
in  the  bibliography.  This  theory  gives  a  list  of  hypotheses  from  which  one  may  prove  that 
optimal  solutions  exists,  that  Lagrange  multipliers  may  be  used  to  enforce  the  constraints, 
and  also  error  estimates  for  approximations  to  optimal  solutions.  We  then  have  shown 
how  problems  from  disparate  areas,  i.e.,  fluids,  nonlinear  elasticity,  and  superconductivity, 
may  be  analyzed  using  our  abstract  theory. 


Another  area  of  research  has  been  in  the  feedback  control  of  unsteady,  viscous,  incom¬ 
pressible  flows.  Details  are  given  below  and  in  papers  [16]  and  [18]  listed  in  the  bibliography. 
Our  main  result  is  to  show  how  a  simple  feedback  law  ran  be  effectively  used  to  reduce 
the  size  of  the  oscillations  in  the  lift  coefficient  in  flow  about  a  cylinder  This  problem  is 
of  interest  in  its  own  right,  e.g.,  in  flows  about  rotor  shafts  in  helicopters,  and  is  also  a 
prototype  for  other,  more  complex  problems. 
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We  have  also  engaged  in  a  study  of  models  and  algorithms  for  macroscopic  supercon -  "od«ie 


.  Vcr 


2 

Diet 

-v-s 

.  c-.l 

ductivity.  Details  are  given  in  papers  [7]  and  [11]  listed  in  the  bibliography.  One  major 
accomplishment  in  this  regard  is  the  development,  analyses,  and  computation  of  a  model 
for  superconducting  thin  films  having  variable  thickness(7].  We  have  also  computed  using 
this  model.  Two  important  features  of  the  model  is  that  it  is  much  simpler,  and  thus 
cheaper  to  compute  with,  that  a  full- hedged  3D  Ginzburg-Landau  model  and  yet  it  ran 
account  for  vortex  pinning  phenomena  due  to  narow  regions  in  the  film. 


Least-squares  finite  element  methods  for  incompressible,  viscous  flows 

Least-squares  finite  element  methods  for  the  numerical  approximation  of  solutions  of 
partial  differential  equations  have  recently  being  gaining  much  attention  in  the  engineering 
and  mathematical  communities.  In  particular,  for  incompressible  viscous  flow  simulations, 
such  methods  applied  to  first-order  formulations  are  gaining  ever  increasing  interest  and 
acceptance.  Some  of  the  reasons  that  these  methods  are  attractive  in  the  Navier-Stokes 
setting  are: 

•  the  choice  of  approximating  spaces  is  not  subject  to  the  LBB  condition  that  arises  in 
Galerkin  mixed  finite  element  approximations; 

•  a  single  approximating  space  can  be  used  for  all  variables; 

•  solution  methods  can  be  devised  that  require  no  matrix  assemblies,  even  at  the  element 
level; 

•  used  in  conjunction  with  an  appropriate  linearization  method,  e.g.,  a  Newton,  re¬ 
sults  in  symmetric,  positive  definite  linear  systems,  at  least  in  the  neighborhood  of  a 
solution; 

•  used  in  conjunction  with  properly  implemented  continuation  (with  respect  to  the 
Reynolds  number)  techniques,  a  solution  method  can  be  devised  that  will  only  en¬ 
counter  symmetric,  positive  definite  linear  systems; 

•  standard  and  robust  iterative  methods  for  symmetric,  positive  definite  linear  systems 
can  be  used; 

•  no  artificial  boundary  conditions  for  the  vorticity  need  be  introduced  at  boundaries 
at  which  the  velocity  is  specified;  and 

•  accurate  vorticity  approximations  are  obtained. 

The  goals  of  our  research  have  been  to  develop,  implement,  and  analyze  least-squares 
finite  element  methods  for  the  Stokes  and  Navier-Stokes  equations.  In  particular,  we 
have  developed  practical  methods  that  can  rigorously  be  shown  to  be  optimally  accurate. 
As  a  result,  least-squares  finite  element  methods  are  now  demonstrably  superior  to  other 
methods  for  incompressible  flow  simulations  in  that  they  require  less  storage,  can  be  solved 
for  more  efficiently,  and  for  the  same  cost  of  simulation,  yield  more  accurate  results.  Here, 
we  will  describe  some  of  our  results  in  conjunction  with  the  Stokes  and  Navier-Stokes 
equations.  Details  may  be  found  in  [l]-[5]. 
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The  best  least-squares  finite  element  method  for  incompressible  flow  calculations  is 
based  on  the  velocity-vorticity-total  pressure  formulation.  For  the  generalized  linear  2D 
Stokes  equations  (we  will  discuss  the  3D  and  nonlinear  cases  below),  we  have  the  system 
of  first-order  differential  equations 

curia/  +  grad  p  =  fi  ] 

curl  n  —  u  =  fi  >  in  Q 
div  u  =  fz  J 

for  the  vorticity  u,  the  velocity  u,  and  the  total  pressure  p.  (The  pressure  itself  is  easily 
determined  from  the  total  pressure  and  the  velocity.)  In  the  fluids  setting,  the  forcing 
functions  /a  and  fo  vanish.  We  consider  (and  contrast)  two  types  of  boundary  conditions: 

u  *  U  on  T  (BCl) 


or 

u  •  n  =  Un.  and  p  =  P  on  T .  (BC2) 

We  assume  that  the  data  satisfy  all  necessary  compatibility  conditions,  e.g.,  for  the  bound¬ 
ary  condition  (BCl), 

f  hdTl=  [  U-ndT. 

J  n  Jr 

For  simplicity,  we  will  consider  homogeneous  boundary  conditions,  i.e.,  U  =  0  or  Un  **  0 
and  P  ss  0. 

In  order  to  facilitate  our  discussion,  we  introduce  some  Sobolev  space  notation.  First, 
for  a  positive  integer  m,  we  have  the  space 

Hm(Q)  ={  set  of  functions  such  that  all  partial  derivatives 

of  order  <  m  are  square  integrable  } 


and  the  associated  norm  for  a  function  g  belonging  to 

K _h2 _  Y-  r  fd(mi+mi)g\2 


2  _  V  f  (#m™9V  i0 


m  x  £  m 


For  example 


Ul-f  srd fi, 

Ja 

r  2  2  T 

+  (f?)  +SS  \dn  =  \\dgldx\\l  +  <lidgiay\\l  +  Ml,  “d 
llfllli  =  W&s/dx'Wl  +  nava*0»n3  +  II^Wllo  +  llsll?  • 
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Our  analyses  of  the  Stokes  problem  uses  the  Agmon-Douglis-Nirenberg  (ADN)  theory 
for  elliptic  systems  of  partial  differential  equations.  Of  particular  interest  to  use  are  the 
following  observations  about  the  ADN  theory. 

•  The  ADN  theory  yields  a  priori  estimates  for  solutions  of  elliptic  boundary  value 
problems. 

•  The  norms  appearing  in  the  estimates  are  chosen  so  that  the  differential  operator 
and  boundary  condition  operator  satisfy  a  certain  precise  condition  known  as  the 
complementing  condition. 

•  If  the  differential  operator  is  elliptic,  and  the  complementing  condition  is  satisfied, 
we  call  the  system  of  partial  differential  equations  and  boundary  conditions  an  ADN 
system. 

•  For  the  same  system  of  differential  equations,  different  boundary  conditions  may  re¬ 
sult  in  the  usage  of  different  norms  within  the  ADN  theory,  i.e.,  a  system  of  partial 
differential  equations  and  boundary  conditions  may  be  an  ADN  system  with  respect 
to  different  norms  than  the  same  partial  differential  equations  with  different  boundary 
conditions. 

•  The  correct  ADN  norms  are  related  to  the  principal  part  of  the  differential  operator, 
e.g.,  the  principal  part  of  the  operator  determines  the  well-posedeness  of  the  problem. 

For  the  pressure-normal  velocity  boundary  condition  (BC2),  the  principal  part  of  the 
Stokes  operator  is  given  by 

curlu/  +  gradp ,  curlu,  and  divu. 

Only  first  derivatives  appear  in  the  principal  part.  As  a  result,  one  has  that  ail  variables 
have  the  same  differentiability  properties.  The  principal  part  operator  along  with  the 
boundary  conditions  uncouple  into  two  well-posed  problems: 

curlu/  -r  grad p  =  f;  in  H 
p  —  P  on  r 


and 


curl  u  =  f 2  and  div  u  =  fz  in  Q 
u  •  n  =  Un  on  f . 


The  ADN  a  priori  estimate  relevant  to  least-squares  methods  for  the  pressure-normal 
velocity  boundary  condition  is  given  by 


Mi  +  Iblli  +  Mi  <  C7(lfi||o  +  II/2II0  +  II/3II0) . 

Note  that  all  norms  on  the  solution  are  the  same,  and  that  also  all  the  data  is  measured 
in  just  one  norm. 
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If,  for  the  velocity  boundary  condition  (BCl),  one  chooses  the  same  principal  part 
as  that  chosen  above  for-  (BC2),  we  see  that  the  principal  part  operator  along  with  the 
boundary  condition  uncouple  into  the  two  problems 

curl  u;  +  grad  p  =  f i  in  0 


and 


curl  u  =  f 2  and  div  u  =  fz  in  fl 

u  =  U  on  T . 


These  problems  are  not  well- posed,  i.e.,  the  first  is  under  determined  (not  enough  boundary 
conditions),  the  second  is  overdetermined  (too  many  boundary  conditions).  The  ADN  the¬ 
ory  tells  us  that  the  principal  part  of  the  Stokes  operator  with  velocity  boundary  conditions 
is  given  by 

curlur  +  gradp,  -w  +  curlu,  and  divu 

so  that  the  principal  part  is  the  whole  Stokes  operator.  The  second  operator  in  the  principal 
part  implies  that  u  and  u  cannot  have  the  same  differentiability  properties.  The  ADN  a 
priori  estimate  relevant  to  least-squares  methods  for  the  velocity  boundary  condition  is 
given  by 

IMIi  +  Iblli  +  Mi  <  cm  Bo  +  0/alli  +  ll/sili)  - 

Note  that  different  components  of  the  solution  are  measured  in  different  norms  and  that 
the  different  components  of  the  data  are  also  measured  in  different  norms. 

Note  the  consistency  achieved  by  the  ADN  theory.  If  u  has  two  square  integrable 
dervatives,  then  u  and  p  have  one  square  integrable  derivative.  Then  the  combination 
curl cj -f- grad p,  i.e.,  fi,  should  be  merely  square  integrable,  and  the  combinations  curl  u-u/ 
and  div  u,  i.e.,  fi  and  fi,  respectively,  should  have  one  square  integrable  derivative.  These 
are  exactly  the  norms  appearing  in  the  a  priori  estimate. 


If  one  uses  the  same  norm  for  all  unknowns  (and  also  the  same  norm  for  all  the  data), 
then,  in  the  velocity  boundary  condition  case,  the  Stokes  system  is  not  an  ADN  system, , 
i.e.,  the  system  is  not  well-posed  with  respect  to  those  norms. 

A  least-squares  functional  can  be  set  up  by  summing  up  the  squares  of  the  residuals 
of  the  equations: 

J(vl,p,oj)  =  ||curlw  +  gradp-fi||2  +  ||curiu-u/  — /2||2  +  ||divu-/3[|2. 


Naturally,  one  asks  the  question:  what' norms  should  be  used  to  measure  the  size  of  the 
residuals?  An  answer: 

o  if  one  uses  the  norms  indicated  by  the  ADN  theory, 
o  and  if  one  also  uses  a  conforming  finite  element  method, 

o  then,  from  a  practiced  point  of  view ,  optimally  accurate  solutions  are  obtained  for  all 
variables 
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o  and,  from  a  mathematical  point  of  view,  the  analysis  of  errors,  e.g.,  the  derivation  of 
rigorous  error  estimates,  is  completely  straightforward. 


For  the  normal  velocity- pressure  boundary  condition  case  (BC2),  the  ADN  theory 
suggests  the  use  of  the  least-squares  functional 

-  ||curlu/4-gradp- +||curlu-w-  /2||J  +  ||divu-/3||o 

=  f  |curlu/  -i-  gradp  -  fi|2  dd  -r  [  (curl  u  -  u  -  ft)2  dCl  +  f  (div  u  -  /3)2  dfl . 

Jo  Jo  Jo. 

Note  that  this  functional  involves  at  most  products  of  first  derivatives.  The  least-squares 
principle  is  then  given  by: 

seek  (u,p,oj)  such  that  J  is  minimised  over  an  appropriate  class  of  functions  V. 

The  function  class  V  consists  of  HX(Q)  velocity,  pressure,  and  vorticity  fields,  constrained 
by  boundary  conditions,  etc.  The  Euler-Lagrange  equation  for  the  least  squares  principle 
are  given  by: 

seek  (u,  p,  a/)  6  V  such  that  5((u,  p,  w) ,  ( v,  q,  0)  =  ^(v,  g,  0  for  all  ( v,  q,  f )  €  V , 
where 

B((n,p,u) ,  (v,g,£))  =  /  (curl u  4- grad p)  •  (curl £  -f  grad q)  dVt 

Jo 

+  [  (curl  u  —  uj)  (curl  v  —  £)  dP.  +  /  (div  u)  (div  v)  dCl 

Jo  Jo 

and 

F(v, <hO~  /  fi  •  (curl  f  +  grad  gjdfi-f-  f  /2(curiv-f)dfl  +  /  /jdivvdfi. 

Jo  Jo  Jo 

Conforming  finite  element  approximations  are  defined  in  the  usual  manner.  One 
chooses  a  conforming  finite  element  approximating  space  V\  i.e.,  the  finite  element  func¬ 
tions  for  all  variables  have  one  square  integrable  derivative.  Then,  one  poses  the  problem: 

seek  (u h,pf%,uh)  6  Vh  such  that 

5((u\p\w*),  (vh,qh,£h))  =^'(v/*,g,^)  for  all  (v\g\f*)  €  V*. 

This  problems  is  equivalent  to  a  linear  algebraic  system  having  a  symmetric,  positive 
definite  coefficient  matrix.  Standard  finite  element  methodology,  i.e.,  based  on  the  Lax- 
Milgram  theorem,  can  be  used  to  derive  optimal  error  estimates.  For  example,  if  piecewise 
linear  polynomials  are  used  for  all  variables  (and  the  exact  solution  is  sufficiently  smooth), 
one  finds  that 

!|u  -  u%  +  ||p  -  ph Hi  +  Ik  -  o//l||1  =  0(h) 

II u  -  uH||o  +  ||p  -  p*||o  +  Ik  -  uh\\o  =  0(h ?) ; 
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if  piecewise  quadratic  polynomials  are  used,  one  finds  that 


ll'u  -  u‘n,  +  lb  -  A  +  ||u  -  =  O(h') 

l|u  -  u*||o  +  lb  -  p'ilo  +  ||w  -  <J*8o  =  0(h3) . 

Note  that  all  variables  are  approximated  by  the  same  finite  element  functions,  all  variables 
are  optimally  approximated,  and  conforming  finite  element  approximations  for  all  variables 
are  required  to  be  merely  continuous  across  element  edges. 


For  the  velocity  boundary  condition  case  (BCl),  the  ADN  theory  suggests  the  use  of 
the  least-squares  functional 

AC(u ,p,w)  =  ||curla/-i-gradp-  fij||  +  ||curlu-u/-/2||f  -f  |{div  u  —  /3i|f 
=s  /  (curio;  -f-  gradp  —  fi|2  dCl  +  /  |grad  (curlu  —  cj  —  /2)|2 

Jcx  Jci 

■f  J  (curlu -w  -  /2)2  dCl  +  j  (|grad  (divu  -  fz)\2  +  (divu  -  /2)2j  dfl 


Note  that  this  functional  involves  products  of  second  derivatives  of  the  velocity.  The  least- 
squares  principle  is  then  given  by: 


seek  (u,  p,  o ;)  such  that  AC  is  minimized  over  an  appropriate  class  of  functions  W. 

The  function  class  W  now  consists  of  Hl(Q)  pressure  and  vorticity  fields  and  H~{ Q) 
velocity  fields,  constrained  by  boundary  conditions.  The  Euler-Lagrange  equation  for  this 
least  squares  principle  again  has  the  form 

seek  (u,  p,  u)  6  W  such  that  5((u,p,w),  (v,g,f))  =  f(vtq,£)  for  all  (v,g,f)  €  W, 

where  now  B  involves  products  of  second  derivatives  of  u  and  v. 

Conforming  finite  element  approximations  are  again  defined  in  the  usual  manner.  One 
chooses  a  conforming  finite  element  approximating  space  Wh,  i.e.,  one  chooses  the  finite 
element  functions  for  approximating  the  pressure  and  vorticity  so  that  their  derivatives 
are  square  integrable  and  finite  element  functions  for  approximating  the  velocity  so  that 
their  second  derivatives  axe  square  integrable.  Then,  one  poses  the  approximate  problem: 

seek  ( n,l,ph,uh )  6  W*  such  that 

5((u\p\wA),  (v*,  $*,£*))  =  .F(v\g,fA)  for  all  (vh,qh,€h)  €  WA. 

This  problems  is  again  equivalent  to  a  linear  algebraic  system  having  a  symmetric,  positive 
definite  coefficient  matrix.  Standard  finite  element  methodology,  i.e.,  based  on  the  Lax- 
Milgram  theorem,  can  be  used  to  derive  optimal  error  estimates  whenever  conforming 
finite  element  spaces  are  used.  However,  the  method  is  not  practical.  The  requirement  that 
finite  element  velocity  approximations  possess  two  square  integrable  derivatives  forces  one 
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to  use  finite  element  functions  that  are  continuously  differentiable  across  element  edges. 
(Incidentally,  if  one  is  willing  to  use  continuously  differentiable  velocity  approximations, 
one  might  as  well  do  least-squares  on  the  primitive  variable  formulation!) 

Thus  we  are  faced  with  the  following  dilemma.  If,  for  velocity  boundary  conditions, 
we  use  a  least-squares  functional  based  on  the  ADN  norms  we  are  led  to  a  computational 
method  requiring  continuously  differentiable  velocity  approximations,  i.e.,  we  get  an  im¬ 
practical  method.  If  instead  we  use  the  more  practical  functional  J  that  works  easily 
amd  optimally  for  the  normsd  velocity/pressure  boundary  conditions,  we  get  non-optimal 
approximations.  Naturally,  one  may  ask  the  following  question: 

is  there  a  way  to  use  the  simpler  and  more  practical  norms  of  the  functional  J  and 
still  get  optimally  accurate  approximations? 

An  affirmative  answer  is  found  by  introducing  mesh-dependent  weights  in  the  least-squares 
functional.  The  residual  norms  we  would  like  to  use  are  given  by 

j] curl w  +  grad p  -  fi||o  ,  ||curlu  —  u  -  /2||o  ,  and  ||divu  -  /3||0 . 

The  residual  norms  the  ADN  theory  would  like  us  to  use 

|| curl w  +  grad p  -  fx||o  ,  ||curlu -u;  -  /2||x ,  and  ||divu  - /3||i . 

A  standard  inverse  inequality  for  finite  element  functions  is  given  by 

n«*iii  <  ca-viio  , 

where  h  is  an  appropriate  measure  of  the  grid  size.  This  suggests  that,  for  finite  element 
functions,  one  can  “simulate”  the  norm  ||qA||i  by  n-1||gA|jo .  Thus,  the  weighted  least 
squares  functional  is  given  by 

s7/i(u,p,w)  =  llcurlu/  +  gradp  -  fi||o  +  /i“2||curl  u  -  w  -  /2||g  +  A~2||divu  -  /3||o 

=  /  [curio;  +  gradp  —  fi|2  dfl 

Jn 

+  jp  J  (curl  u  —  (J  -  f2)2  d&  +  ~  J  (div  u  -  /3)2  dfi . 

This  functional  uses  the  norms  that  lead  to  a  simple,  easy  to  implement  algorithm  for  which 
one  can  use  merely  continuous  finite  element  functions  for  all  variables  and  for  which  one 
still  obtains  a  symmetric,  positive  definite  discrete  linear  system.  This  functional  also  leads 
to  optimally  accurate  approximations.  Note  that  optimal  accuracy  is  achieved  with  respect 
to  norms  related  to  the  ADN  theory  for  velocity  boundary  conditions. 

~  The  finite  element  analyses  indicate  that  one  may  use  polynomials  of  one  degree  lower 
for  the  pressure  and  vorticity  than  one  uses  for  the  velocity.  If  we  use  continuous  piecewise 
quadratic  polynomials  for  the  velocity  approximations  and  piecewise  linear  polynomials 
for  the  pressure  and  vorticity  approximations,  we  get  the  estimate 

llu  - u*n, +  ||p-/||o  +  II*  -  *ic  =  o (h*) . 
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This  estimate  is  optimal  with  respect  to  the  finite  element  functions  used.  The  theory  also 
says  that  one  may  use  the  same  degree  polynomials  for  all  variables.  For  example,  if  one 
uses  continuous  piecewise  quadratic  polynomials  for  all  variables,  one  again  obtains  the 
above  error  estimate.  In  this  case,  the  above  estimate  is  not  optimal  for  the  pressure  and 
vorticity. 

In  three  dimensions,  the  velocity- vorticity-pressure  formulation  of  the  Stokes  equations 
with  velocity  boundary  conditions  is  given  by 

curio;  +  gradp  =  fi 
curl  u  —  ijJ  =  f2 
div  u  =  /3 


and 


u  =  U  on  T , 


where  now  the  vorticity  a;  is  a  vector  valued  function.  We  immediately  have  a  problem:  we 
have  seven  equations  in  seven  unknowns  so  that  the  system  cannot  be  elliptic.  We  avoid 
this  problem  by  adding  an  extra  variable  6  and  a  seemingly  redundant  equation  to  get  the 
elliptic  system  of  8  equations  in  8  unknowns: 


curio;  -i-  grad  p  =  fi 

divo;  =  —  divf2 
curl  u  -f-  grad  <p  -  u  =  f2 
div  u  =  fz 


in  Cl . 


One  can  easily  show  that  <6  =  0.  Analyses  yield  the  same  results  as  in  the  2D  case.  Algo¬ 
rithmically,  one  can  completely  ignore  p,  but  one  cannot  ignore  the  redundant  equation. 
Thus,  we  can  discretize  the  problem 


curio;  -f-  gradp  =  f; 

divo;  = —div  fj  .  „ 

curl  u  —  u  —  f2 
divu  =  /3 

u  =  U  on  T 

using  the  weighted  least  squares  functional 

Jh(u,p,u)  =  |j curio;  4-  gradp  -  fi||J  +  || divo;  4-  div  f2||| 

+  h“2||curl  u  -  u  -  f2||o  4-  h_2||divu  -  /3||s 

=  /  |curlo;  4-  gradp  -  fi|2  dCl  4-  /  (divo;  4-  divf2)2  dCl 

Jn  J  n 

4-  i  J  (curl  u  -  o)  -  f2)2  dCl  +  ^p  j  (div  u  -  /3)2  dCl 
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The  above  discussion  can  be  extended  to  the  nonlinear  Navier-Stokes  equations 


i/curlo;  -f-  gradp  4-  a/  x  u  =  fi 

curl  u  —  u  =  0  ► 

div  u  =  0 


in  Q 


and 


u  =  U  on  T . 


The  weighted  least  squares  functional  in  this  case  is  given  by 


Jh.(u,p,u )  ss  v~-\\i/cur\u  4-  gradp  +  £J  x  u  -  fxj(o  4-  i/“2||divu;||Q 
+  /i-2||curlu-u/||o  4-  /i”2||divu||o 

=  f  j^curlw  -f  gradp  4-  x  u  —  fx|2  dTt-r  X;  f  jdiv tu[2  dfi 
y‘  J n  u~  Ju 

+  Z*  [  (curl  u- a;)2  dfl  +  f  (divu)2  dCl. 

J fl  Jn 


All  theoretical  results  for  the  Stokes  equations  can  also  be  derived  in  the  context  of  the 
nonlinear  Navier-Stokes  equations.  Algorithmically,  one  must  choose  a  method  for  lineariz¬ 
ing  the  equation,  if  one  uses  Newton’s  method,  then  in  the  neighborhood  of  a  solution, 
the  Hessian  matrix  is  not  only  symmetric,  but  is  also  positive  definite. 


The  analysis  of  the  nonlinear  case  required  and  extension  of  the  Brezzi-Rappaz-Raviart 
theory  for  the  approximation  of  nonlinear  problems.  That  theory  addresses  problems  of 
the  following  type.  We  let  X  and  Y  denote  Banach  spaces,  A  a  compact  interval  of  the 
real  line,  T  :  Y  —*  X  a  linear  operator,  and  G  :  X  x  A  — ►  Y  a  nonlinear  operator.  Then, 
for  A  €  A,  we  seek  4>  £  X  such  that 


ip+TG(j;,  A)  =0. 

Approximations  are  defined  as  follows.  Let  Xh  C  X  and  let  Th  :  Y  Xh  be  a  linear 
operator.  Then,  for  A  €  A,  we  seek  iih  €  XK  such  that 

1>h  +  ThG{1>h,  A)  =0. 

Under  certain  conditions  on  the  G,  it  turns  out  that  estimates  for  the  error  || $  —  i>h\\x 
can  be  determined  from  estimates  for  j|X  —  Tfc||.  This  theory  has  been  widely  applied; 
for  example,  it  can  be  used  to  determine  error  estimates  for  Galerkin  mixed  finite  element 
approximations  of  the  primitive  variable  formulation  of  the  Navier-Stokes  equations.  It 
can  also  be  applied  directly  to  conforming  least  squares  finite  element  approximations  of 
the  velocity-vorticity-pressure  formulation  of  the  Navier-Stokes  equations  based  on  the  use 
of  the  norms  suggested  by  the  ADN  theory.  In  these  cases,  the  operator  T  is  merely  the 
solution  operator  of  the  corresponding  Stokes  problem  and  an  estimate  for  \\T  —  T*||  can 
be  obtained  from  the  analyses  of  the  analogous  least-squares  finite  element  approximation 
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of  the  Stokes  equations.  One  finds  that  the  errors  in  approximations  of  solutions  of  the 
nonlinear  Navier-Stokes  equations  behave  in  the  same  manner  (with  respect  to  the  grid 
size)  as  do  the  errors  in  the  analogous  approximations  of  solutions  of  the  linear  Stokes 
equations. 

Least-squares  finite  element  methods  for  the  Navier-Stokes  equations  based  on  the 
use  of  weighted  functionals  do  not  fit  into  the  Brezzi-Rappaz-Raviart  framework.  The 
problem  fits  into  the  form  i)  4-  TG(i),  A)  =  0  using  the  space  X  suggested  by  the  ADN 
theory.  However,  although  the  dicsrete  problem  is  of  the  form  +  ThG(i>h,  A)  =  0, 
the  discrete  space  XK  £  X,  i.e.,  we  are  within  the  realm  of  nonconforming  finite  element 
methods.  We  have  extended  the  Brezzi-Rappaz-Raviart  theory  to  the  nonconforming  case 
and  have  shown,  for  least-squares  finite  element  methods  for  the  Navier-Stokes  equations 
based  on  wieghted  functionals  that  optimal  error  estimates  are  obtained. 

We  give  the  results  of  some  computational  experiments  merely  to  show  the  necessity 
of  introducing  the  weights  into  the  functional.  For  these  figures,  we  use  piecewise  linear 
pressures  and  vorticities  and  piecewise  quadratic  velocities.  (Similar  results  can  be  ob¬ 
tained  for  piecewise  quadratic  velocities  and  pressures.)  In  the  figures  below,  there  are  two 
graphs  in  each  frame.  The  one  with  the  steeper  slope  corresponds  to  the  use  of  the  weighted 
functional,  while  the  other  one  corresponds  to  the  use  of  the  unweighted  functional.  Thus, 
we  see  that  the  use  of  the  unweighted  functional  results  in  a  loss  of  accuracy.  Morever,  one 
can  verify  that  the  slope  of  the  graphs  correponding  to  the  use  of  the  weighted  functional 
indicate  that  those  approximations  converge  at  an  optimal  rate. 


Hi  Srrar:  u- velocity  SI  Srrsr:  r-v« lacisy 


Analysis  of  a  shape  control  problem  for  the  Navier-Stokes  equations 

Shape  controls  are  of  great  interest  in  the  control  of  fluid  flows.  For  example,  the 
problem  of  the  optimal  design  of  wings  is  a  shape  optimal  control  problem.  There  have 
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been  substantial  analyses  of  linear  shape  control  problems  for  partial  differential  equations. 
However,  for  the  nonlinear  Navier-Stokes  equations  of  incompressible  flow,  no  rigorous 
analyses  have  been  previously  carried  out.  We  have  given  the  first  such  analyses  in  the 
context  of  a  two-dimesional  flow  in  a  channel-like  domain.  Our  rigorous  analyses  include 
showing,  in  the  context  of  a  drag  minimization  problem,  that  optimal  solutions  exist, 
showing  that  the  Lagrange  multplier  may  be  used  to  enforce  constraints,  deriving  an 
otpimality  system  from  which  optimal  states  and  co-states  may  be  deduced,  deriving  the 
shape  gradient,  and  finally,  deriving  error  estimates  for  finite  element  approximations  of 
optimal  states.  Details  will  be  provided  in  (14|  and  [15 j. 

Our  prototype  problem  is  that  of  two-dimensional  incompressible  flow  of  a  viscous 
fluid  passing  through  a  channel  having  a  bump  on  the  lower  wall;  see  the  figure  below. 
(Our  results  can  be  extended  to  a  variety  of  other  problems.)  We  describe  the  bump 
through  the  function  y  =  ct(x);  note  that  the  extent  of  the  bump  [xi,  xr\  is  fixed  and  that 
a(xi)  =  a(xr)  =  0  are  imposed  as  constraints  on  ct(z).  We  also  require  that  cc(x)  <  L. 
Thus,  the  bump  is  described  by 

Ta  =  { (*,  y )  €  [zi,  xr]  x  [0,  L\  :  y  =  a(x)  }  . 

We  let  T  denote  the  boundary  of  the  flow  domain  Oa;  of  course,  ra  C  I\  The  admissible 
family  of  curves  Fa  is  defined  by 

Had.  =  {d  €  C0'l{xi,xr)  :  0<a(x)<L,  |a'(r)|  </3Vx  ~  [xi,xr\ ,  a(x*)  =  a{xr)  =  0}  , 

where  the  constant  ,3  >  0  is  chosen  so  that  Uad  7=  0  and  C°’l(xi,xr)  denotes  the  Lipschitz 
continuous  functions  denned  on  [xt,xr ]. 


v=L 


The  state  equations  are  defined  on  P.a  for  each  a  €  Uad-  They  are  given  by  the 
Navier-Stokes  system: 


and 


-i/Au  +  u  •  V  U--r  Vp  =  f  in  , 
V  •  u  =  0  in  , 


fg  onT /Ta 
\0  on  r a  , 
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where  f  and  g  ar  given  fucntions.  (Other  boundary  conditions  on  r/ra  can  also  be  treated.) 
The  objective  functional  we  treat  is  the  dissipation  functional 

J(a)  =  J(Cla,  ua)  =  v  f  Vua  :  Vua  dx , 

where  u(a)  is  a  solution  of  the  state  equations  in  the  domain  Cla.  (We  could  also  consider 
other  functionals  such  as  those  involving  the  matching  of  the  flow  to  some  given  flow.) 

The  extremal  problem  we  consider  is  then  given  by: 

min  J{ct)  subject  to  (ua,pa)  satisfying  the  Navier-Stokes  system  over  Qa . 

ot€Umj 

Our  first  main  result  is  to  show  that  this  extremal  problem  has  a  solution.  We  then  show 
that  the  use  of  the  Lagrange  multiplier  rule  to  enforce  the  contraints,.  i.e.,  the  Navier- 
Stokes  system,  is  rigorously  justified.  Subsequently,  we  derive  the  necessary  conditions 
that  optimal  controls  and  states  must  satisfy.  These  are  found  by  requiring  that  variations 
in  an  appropriate  Lagrangian  vanish.  Variations  with  respect  to  the  Lagrange  multipliers 
yields,  of  course,  the  state  constraint  equations,  i.e.,  the  Navier-Stokes  system.  Variations 
in  the  state  (u ,p)  yield  the  co-state  equations 

-t/Aq  —  u  •  Vq  +  q  •  (Vq)r  +  Vr  =  -2i/Au  in  Qa  , 

V  •  q  ss  0  in  flQ , 

and 

q  ss:  0  on  T . 

Note  that,  as  usual,  these  are  the  adjoint  of  the  linearized  Navier-Stokes  system  where 
the  linearization  is  about  the  state  u.  Variations  in  the  control  function  a  can  be  used  to 
determine  the  shape  gradient  of  the  functional  J,  i.e.,  roughly  speaking,  the  derivative  of 
J  with  respect  to  a.  This  is  given  by 

gradj  =  i/Vu  :  Vu  —  i/V u  :  Vq  +  pV  •  q  —  (n  •  Vu)  •  (2u(n  •  Vu)  —  u(n  •  Vq)  4-  raj  . 

Thus,  we  see  that,  for  a  given  function  a,  the  shape  gradient  gradj  depends  on  the  state 
(u,p)  and  the  adjoint  state  (q,  r)  determined  from  the  state  and  co-state  systems  posed 
over  the  domain  Qa,  i.e.,  for  a  given  a,  gradj(a )  =  gradj(u.a,pa,  qa,ra).  (The  above 
presentation  is  formal;  the  rigorous  derivations  we  have  done  are  with  respect  to  weak 
formulations  and  solutions  posed  over  appropriate  function  spaces.) 

A  typical  optimization  algorithm  would  use  the  shape  greadient  gradj ,  as  well  as  J 
itself,  to  update  a  guess  for  the  shape  function  a.  As  we  have  seen,  we  may  compute  the 
shape  gradient  by  for  any  iterate  a  by  solving  the  state  equations,  i.e,  the  Navier-Stokes 
system,  and  the  co-state  equations  with  respect  to  the  domain  fta.  Then,  these  solutions 
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axe  used  to  evaluate  the  functional  and  the  shape  gradient  of  the  functional  J  using  the 
above  formulas.  Thus,  the  above  equations  completely  determine  the  information  needed 
by  any  optimization  algorithm  requiring  gradient  information.  Our  main  contributions  in 
this  respect  are  to  derive  the  sjhape  gradient  formula  and  to  rigorously  justify  all  the  steps 
performed  in  that  derivation. 

Of  course,  in  general  one  can  only  solve  the  state  and  co-state  systems  in  an  approxi¬ 
mate  manner.  To  this  end  we  have  defined  finite  element  methods  for  the  approximation 
of  the  state  and  co-state  equations  (for  a  given  value  of  a),  and  derived  optimal  error 
estimates  for  the  finite  element  approximations.  We  have  also  defined  approximations  for 
the  stress  and  adjoint  stress  vectors  (these  appear  in  the  definition  of  the  shape  gradient) 
and  derived  optimal  error  estimates  for  these  quantities  in  appropriate  norms.  Thus,  we 
also  have  error  estimates  for  approximations  of  the  shape  gradient. 


Finite  dimensional  approximation  of  a  class  of  nonlinear  optimal  control  prob¬ 
lems 

The  need  to  solve  optimization  and  control  problems  arises  in  many  settings.  Although 
in  some  cases  these  problems  may  be  easily  solved,  either  analytically  or  computationally,  in 
many  other  cases  substantial  difficulties  are  encountered.  For  example,  candidate  optimal 
states  and  controls  may  belong  to  infinite  dimensional  function  spaces  and  one  may  have 
to  minimize  a  nonlinear  functional  of  the  state  and  control  variables  subject  to  nonlinear 
constraints  that  take  the  form  of  a  system  of  partial  differential  equations  whose  solutions 
are  in  general  not  unique.  Our  goal,  which  we  have  reached  (see  [13]  for  details),  was 
to  construct,  analyze,  and  apply  a  framework  for  the  approximate  solution  of  many  such 
problems.  The  setting  for  our  framework  is  a  class  of  nonlinear  control  or  optimization 
problems  which  is  general  enough  to  be  of  use  in  numerous  applications.  The  major  steps 
in  the  development  and  analysis  of  our  framework  axe  as  follows: 

•  define  an  abstract  class  of  nonlinear  control  or  optimization  problems; 

•  show  that,  under  certain  assumptions,  optimal  solutions  exist; 

•  show  that,  under  certain  additional  assumptions,  Lagrange  multipliers  exist 
that  may  be  used  to  enforce  the  constraints; 

•  use  the  Lagrange  multiplier  technique  to  derive  an  optimality  system  from 
which  optimal  states  and  controls  may  be  deduced; 

•  define  algorithms  for  the  finite  dimensional  approximation  of  optimal  states 

and  controls;  and  ' 

•  derive  estimates  for  the  error  in  the  approximations  to  the  optimal  states 
and  controls. 

Two  of  the  key  ingredients  used  to  carry  out  the  above  plan  axe  the  Tikhomorov  version  of 
the  Lagrange  multiplier  theory  and  the  Brezzi-Rappaz-Raviart  theory  for  the  approxima- 
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tion  of  a  class  of  nonlinear  problems.  In  both  of  these  theories,  certain  properties  of  com¬ 
pact  operators  on  Banach  spaces  play  a  central  role.  We  point  out  that  the  nonuniqueness 
of  solutions  of  the  nonlinear  constraint  equations  deems  it  appropriate  to  employ  Lagrange 
multiplier  principles. 

After  having  developed  and  analyzed  the  abstract  framework,  we  applied  it  to  some 
specific,  concrete  problems.  In  each  case,  we  used  the  abstract  framework  to  analyze  the 
concrete  problems  by  merely  showing  that  the  latter  fit  into  the  former.  The  particular 
applications  we  considered  are: 

•  control  problems  in  structural  mechanics  having  geometric  nonlinearities  that 
are  governed  by  the  von  Karmin  equations; 

•  control  problems  in  superconductivity  that  are  governed  by  the  Gin  r- 
Landau  equations;  and 

•  control  problems  for  incompressible,  viscous  flows  that  are  governed  by  the 
Navier-Stokes  equations. 

In  considering  these  applications,  we  will  purposely  chose  different  types  of  controls  in 
order  to  illustrate  how  these  could  be  accounted  for  within  the  abstract  framework.  In  all 
three  cases,  approximation  were  effected  through  the  use  of  finite  element  methods. 

We  now  give  a  precise  definition  of  the  abstract  class  of  nonlinear  control  or  optimiza¬ 
tion  problems  that  we  have  studied  and  for  which  the  above  results  have  been  obtained. 
We  introduce  the  spaces  and  control  set  as  follows.  Let  G,  X,  and  Y  be  reflexive  Banach 
spaces  whose  norms  are  denoted  by  ||  •  || c,  II  *  ||*y»  and  ||  •  ||y,  respectively.  Dual  spaces 
will  be  denoted  by  (•)*.  The  duality  pairing  between  X  and  X"  is  denoted  by  (*,  -)x;  one 
similarly  defines  (•,  -}y  and  {•,  -)c.  The  subscripts  are  often  omitted  whenever  there  is  no 
chance  for  confusion.  Let  0,  the  control  set,  be  a  closed  convex  subset  of  G.  Let  Z  be  a 
subspace  of  Y  with  a  compact  imbedding.  Note  that  the  compactness  of  the  imbedding 
Z  c  Y  plays  an  important  role.  We  assume  that  the  functional  to  be  minimized  takes  the 
form 

J{v ,  z)  =  A  F{y)  +  A  £(z)  V  (v,  z)  €  X  x  © , 

where  T  is  a  functional  on  X,  £  a  functional  on  0,  and  A  is  a  given  parameter  which  is 
assumed  to  belong  to  a  compact  interval  ACR+.  The  constraint  equation  M(v,z )  =  0 
relating  the  state  variable  v  and  the  control  variable  z  is  defined  as  follows.  Let  N  be  a 
differentiable  mapping  from  X  to  Y,  K  a  continuous  linear  operator  from  ©  to  Y,  and  T 
a  continuous  linear  operator  from  Y  to  X.  For  any  A  €  A,  we  define  the  mapping  M  from 
X  x  0  to  X  by 

M(v, z)=v  +  XTN(v)  +  A TK(z)  V  (u, z)  6  X  x  © . 

With  these  definitions  we  now  consider  the  constrained  minimization  problem: 

min  J{v,z)  subject  to  M(v, z)=0. 

(v  ,x)6Xx& 
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Here,  we  are  seeking  a  global  minimizer  with  respect  to  the  set  {(v,  zj  €  X  x  9  :  M(v,  z)  = 
0}.  Although,  under  suitable  hypotheses,  we  have  shown  that  this  problem  has  a  solution, 
in  practice,  one  can  only  characterize  local  minima,  i.e.,  points  (it,  g)  6  X  x  ©  such  that 
for  some  «  >  0 

J{u,g)  <  »7(v,  z)  V(v,z)  6  X  x  9  such  that  M(v,z )  =  0  and  ||tx  —  v||x  <  e- 

Thus,  when  algorithms  for  locating  constrained  minima  of  J  are  considered,  one  must  be 
content  to  find  local  minima. 

After  having  showed  that  optimal  solutions  exist  and  that  one  is  rigorously  justified 
in  using  the  Lagrange  multiplier  rule,  we  introduce  some  simplifications  in  order  to  render 
the  abstract  problem  more  amenable  to  approximation.  The  first  is  to  only  consider  the 
control  set  ©  =  G.  The  second  is  to  only  consider  Frechet  differentiable  functionals  £(•) 
such  that  the  Frechet  derivative  £'(g)  =  E~lg,  where  E  is  an  invertible  linear  operator 
from  Gm  to  G. 

In  order  to  prove  our  results,  we  will  have  to  introduce  certain  hypotheses.  The  first 
set  of  hypotheses  are  invoked  to  prove  the  existence  of  optimal  solutions.  It  is  given  by: 

(Hi)  inf„€*  F{y)  >  -oo ; 

(H2)  there  exist  constants  a,  f3  >  0  such  that  6{z)  >  o||z|j^  V  z  6  © ; 

(H3)  there  exists  a  (v,z)  €  X  x  0  satisfying  M(v,  z)  =  0 ; 

(H4)  if  — *  u  in  X  and  g^  —  g  in  G  where  {{u^\g^)}  C  X  x  ©,  then 

N{u™)  -  N{u)  in  Y  and  K(g^)  -  K{g)  in  Y ; 

(Ho)  J{‘,  ■)  is  weakly  lower  semicontinuous  on  X  x  © ;  and 

(H6)  if  {(u(n),0(n))}  C  XxQ  is  such  that  {^(u^)}  is  a  bounded  set  in  R  and 
M(u(n\  g^)  =  0,  then  {u^}  is  a  bounded  set  in  X. 

The  second  set  of  assumptions  are  used  to  justify  the  use  of  the  Lagrange  multiplier  rule  and 
to  derive  an  optimality  system  from  which  optimal  states  and  controls  may  be  determined. 
The  second  set  is  given  by: 

(H7)  for  each  z  6  9,  v  ►-*  J(v,z)  and  v  >-*  M(v,  z )  are  Frichet  differentiable ; 

(H8)  z  h-*  £(z)  is  convex,  Le., 

£(7*1  -r  (1  -7)*2)  <  7^(*i)  +  (1-7 )£(*z)  Vzi,z2  €  ©,  V7  €  [0,1]; 
and 

(H9)  for  v  €  X,  N’{v)  maps  X  into  Z. 

In  (H9),  N'  denotes  the  Frechet  derivative  of  N.  A  simplified  optimality  system  may  be 
obtained  if  one  invokes  the  additional  assumption: 

(H10)  ©  a*  G,  and  the  mapping  z>~*  £(z)  is  Frichet  differentiable  on  G. 
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Hypotheses  (HT)-(HIO)  allow  us  to  obtain  a  simplified  optimality  system  for  almost  all 
values  of  the  parameter  *A  6  A.  In  many  cases,  it  is  possible  to  show  that  the  same 
optimality  system  holds  for  all  values  of  A.  The  following  two  additional  assumptions 
which  will  only  be  invoked  in  case  (1/A)  is  an  eigenvalue  of  - TN'(u )  each  provides  a 
setting  in  which  this  last  result  is  valid: 

(Hll)  ifvm  6  **  satisfies  (I  +  A  [;V'(u)]T-)v*  =  0  and  KmT*vm  =  0,  then  v“  «  0 ; 
or 

(Hll);  the  mapping  (v,  z)  v  +  A  TN'(u)v  +A TKz  is  onto  from  X  x  G  to  Y. 

In  order  to  make  the  optimality  system  more  amenable  to  approximation  and  computation, 
we  invoke  the  following  additional  assumption: 

(H12)  £'(g)  =  E~lg ,  where  E  is  an  invertible  linear  operator  from  Gm  to  G  and  g  is 

an  optimal  control  for  the  constrained  minimization  problem  (2.4) . 

Assumptions  (Hl)-(H6)  can  be  used  to  establish  that  optimal  solutions  exist.  With 
the  addition  of  (H7)-(H9),  one  can  show  that  the  Lagrange  multiplier  rule  may  be  used 
to  turn  the  constrained  minimization  problem  into  an  unconstrained  one  and  to  derive 
an  optimality  system  from  which  optimal  states  and  controls  may  be  determined.  By 
also  invoking  (H10),  one  of  (Hll)  or  (Hll*),  and  (H12)  a  simplified  optimality  system 
is  achieved.  In  this  case,  an  optimal  state  u  6  X,  an  optimal  control  g  6  G,  and  the 
corresponding  Lagrange  multiplier  £  €  Ym  satisfy  the  optimality  system 


u  + ATiV(u)  4-  XTKg  =  0  in*, 

£  +  AT* [jV'(u)]*£  -  A TmF(y)  =0  in  V * , 


and 


g  -  EK’S  =  0  in  G. 

In  many  applications  we  have  that  *"  =  Y.  Since  these  spaces  are  assumed  to  be  reflexive, 
we  also  have  that  Y’  =  X.  In  this  case,  we  have  that  both  u  and  f  belong  to  X. 


A  finite  dimensional  discretization  of  the  optimality  system  is  defined  as  follows.  First, 
one  chooses  families  of  finite  dimensional  subspaces  Xh  C  X,  (Ym)h  C  Ym,  and  Gn  C  G. 
These  families  are  parameterized  by  a  parameter  h  that  tends  to  zero.  (For  example,  this 
parameter  can  be  chosen  to  be  some  measure  of  the  grid  size  in  a  subdivision  of  Cl  into 
finite  elements.)  Next,  we  define  approximate  operators  Th  :  Y  —>  Xh,  Eh  :  Gm  -*  Gh, 
and  (T*)a  :  **  —  ( Ym)h .  Of  course,  one  views  T\  Eh,  and  (T*)A  as  approximations  to 
the  operators  T,  E>  and  T*,  respectively.  Note  that  (T*)A  is  not  necessarily  the  same  as 
(TA)\  The  former  is  a  discretization  of  an  adjoint  operator  while  the  later  is  the  adjoint 
of  a  discrete  operator.  Once  the  approximating  subspaces  and  operators  have  been  chosen, 
an  approximate  problem,  or  the  discrete  optimality  system,  is  defined  as  follows.  We  seek 
uK  €  Xh ,  gh  €  Gh,  and  6  ( Ym)H  such  that 

uh  +  A  ThN(uh)  +  A  ThKgh  =  0  in  Xh , 

Sh  +  A  (T*)A(iV,(uA)]-^  -  A  (T*)a^(ua)  =0  in  (Ym)h , 
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and 


gh  -  E*'K'm$h  =  0  in  Gk . 

We  make  the  following  hypotheses  concerning  the  approximate  operators  Th,  (T“)k,  and 
EK: 

(H13)  lim  |J(r  -  rA)y||*  =0  V  y  €  F , 

A— »U 

(H14)  lim  ||  (T*  -  (T-)A)t;||y.  =0  V  v  e  Xm , 

A— *0 

and 

(H15)  lim  ||(5 - Eh)s\\c  =  0  VssG*. 

h—+Q 

We  also  need  the  following  additional  hypotheses  on  the  operators  appearing  in  the  defi¬ 
nition  of  the  abstract  problem: 

(H16)  N  €  C3(X ;  Y )  and  F  6  C3(X;  R) ; 

(H17)  N",  iV ,  T" ,  and  F"  are  locally  bounded,  i. e.,  they  map  bounded  sets  to 

bounded  3ets; 

(H18)  forv  €  X,  in  addition  to  (H9),  i.e.,  N'(v)  6  £(X;  Z )  where  Z  «— »t— ►  Y,  we  have 
that  [iV'(v)j“  €  C(Y’)  Z)  where  Z  — ^  Xm ,  that  for  rj  6  Ym,  [N”(v) ]•  •  tj  6 
£(Ym;  Z),  and  that  for  w  €  X,  F\v)  •  w  6  £(X;  Z);  and 

(H19)  K  maps  G  into  Z. 

Here,  (•)"  and  (•)/,/  denote  second  and  third  Frechet  derivatives,  respectively.  Using  hy¬ 
potheses  (H13)-(H19)  we  can  prove  the  following  results.  Let  (u(A),y(A),£(A))  €  X ,  for 
A  €  A,  be  a  branch  of  regular  solutions  of  the  optimality  system.  Then,  there  exists  a 
5  >  0  and  an  ho  >  0  such  that  for  h  <  Aq,  the  discrete  optimality  system  has  a  unique 
solution  (ua(A),  yA(A),  £a(A))  satisfying 

||(U(A),K(A),?(A))  -  (u‘(A),j‘(A),?‘(A))|U  <  S. 

Moreover, 

Um  ||(t<(A),ff(A),f(A))  -  (tt»(A),/(A),{\A))||*  =  0 
uniformly  in  A  €  A  and  there  exists  a  constant  C,  independent  of  h  and  A,  such  that 
lim  H(»(A),j(A),{(A))  -  (B*(A),j'‘(A),«‘(A))||* 

S  CA{||(T*  -  T)(N(u(\))  +  Kg(X))ftx  + ||(£*  -  £) JC*<(A) ||c 
+ 1|  ((I- )*  -  7~)  ([N'(u(>.m  -  r (u(A)))  ||K. } . 
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These  results  may  be  applied  to  the  derivation  of  error  estimates.  They  basically  state 
that  the  error  in  the  approximate  solution  of  the  optimal  control  problem  is  essentially  the 
same  as  that  in  the  approximations  TK,  (T*)A  and  Eh  to  the  linear  operators  T,  T’,  and 
E.  The  latter  errors  can  usually  be  deduced  using  standard  methodologies,  e.g.,  for  linear 
elliptic  problems. 

The  above  framework  and  analyses  have  been  applied  to  some  concrete  problems,  all  of 
which  feature  constraints  on  admissible  states  and  controls  that  take  the  form  of  a  system 
of  nonlinear  partial  differential  equations.  In  each  application,  we  use  a  different  control 
mechanism  so  that  the  discussion  provided  in  [13]  illustrates  the  treatment  of  a  variety 
of  such  mechanisms.  However,  one  could  use  any  of  the  control  mechanisms  discussed 
in  any  of  the  applications  in  any  other  application,  or  in  fact,  use  any  combination  of 
such  mechanisms.  The  first  application  is  to  distributed  controls  for  the  von  Kdrmdn  plate 
equations.  For  this  application  we  will  use  distributed  controls,  i.e.,  control  is  effected 
through  a  source  term  in  the  governing  partial  differential  equations.  Let  ft  be  a  bounded, 
convex  polygonal  domain  in  R3  and  let  T  denote  the  boundary  of  ft.  The  von  Kdrmdn 
equations  for  a  clamped  plate  are  given  by  (after  suitable  rescaling) 


^20i  +  ;r[02>02l  =  0 


in  ft, 


and 


A202  A  [0i, 0a]  =  A  g  in  ft , 


01  = 


302 

dn 


0  on  r. 


Here,  0t  denotes  the  (sccaled)  Airy  stress  function,  02  the  (scaled)  deflection  of  the  plate 
in  the  direction  normal  to  the  plate,  Xg  is  an  external  load  normal  to  the  plate  which 
depends  on  the  loading  parameter  A,  and  9(-)/5n  the  normal  derivative  in  the  direction  of 
the  outer  normal  to  T.  We  define  the  functional 


9 )  =  3 (0i 1 02 1 9)-^J^  ((tfi  ~  i>io)2  +  (V>2  -  02o)2)  d£l  +  ^  J^g2dn. 

We  then  consider  the  following  optimal  control  problem  associated  with  the  von  Kdrmdn 
plate  equations: 


min  j  0  6  y ,  g  €  ©  }  subject  to  (0,  g)  satisfying  the  von  Kdrmdn  equations , 

where  y  and  @  are  suitable  sets.  In  particular,  functions  in  y  have  square  integrable 
second  derivatives.  For  this  problem,  the  abstract  framework  can  be  used  to  show  that 
optimal  solutions  exist  and  that  the  Lagrange  multiplier  rule  may  be  used,  to  derive  an 
optimality  system  from  which  optimal  states  and  controls  may  be  deduced,  and  to  derive 
optimal  error  estimates  for  conforming  finite  dement  approximations  of  solutions  of  the 
optimality  system. 
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The  second  application  is  the  Neumann  boundary  controls  of  a  simplifie  model  for 
superconductivity.  Let  H“be  a  bounded  open  domain  in  R*,  d  =  2  or  3,  and  let  T  be  its 
boundary.  A  simplified  Ginzburg-Landau  model  for  superconductivity  is  given  by  (after 
appropriate  rescaiings) 

-  A0i  +  (|A|2  -  l)0t  -  V  •  (A02)  -  A  •  V02  4-  A  (02  +  0§)  01=0  in  Cl , 

-A02  +  (|A|2  -  1)02  +  V  •  (A0i)  +  A  •  V0i  +  A  (02  +  02 )  02  =  0  in  Cl , 

n  •  (V0!  +  A02)  =  Xg\  on  T , 
and 

n  •  (V02  —  A0j)  =  A^2  onT. 

Here,  0i  and  02  denote  the  real  and  imaginary  parts,  respectively,  of  the  (rescaled) 
complex-valued  order  parameter,  A  is  a  given  real  magnetic  potential,  <71  and  <72  a-re 
related  to  the  normal  component  of  the  current  at  the  boundary,  and  A  >  0  is  a  “current 
loading”  parameter.  Given  a  desired  state  0<j  =  (0io,  02o)»  we  define  for  any  0  =  (0i,  02) 
and  g  =  (<7i,<7c)  the  functional 

*7(0 , g)  =  ^  ((0i  ~  0io)2  +  (02  -  02o)2)  dCl  +  ^  J^gl  +ai)dr. 

We  then  consider  the  following  optimal  control  problem  associated  with  the  Ginzburg- 
Landau  equations  for  superconductivity: 

mia{*7(0,g)  |  0  €  W,  g  6  ©} 

subject  to  (0, g)  satisfying  the  Ginburg-Landau  equations, 

where  W  and  0  are  suitable  chosen  sets;  in  particular,  functions  belonging  to  W  have 
one  square  integrable  derivative.  For  this  problem,  the  abstract  framework  can  be  used  to 
show  that  optimal  solutions  exist  and  that  the  Lagrange  multiplier  rule  may  be  used,  to 
derive  an  optimality  system  from  which  optimal  states  and  controls  may  be  deduced,  and 
to  derive  optimal  error  estimates  for  conforming  finite  element  approximations  of  solutions 
of  the  optimality  system. 

The  third  aplication  is  to  the  Dirichlet  boundary  control  of  the  Navier-Stokes  equations 
of  incompressible,  viscous  flow.  Let  Cl  denote  a  bounded  domain  in  It*,  d  =  2  or  3,  with 
a  boundary  denoted  by  T.  Let  u  and  p  denote  the  velocity  and  pressure  fields  in  Cl.  The 
Navier-Stokes  equations  for  a  viscous,  incompressible  flow  are  given  by  (after  appropriate 
rescaiings) 

— V •  ((Vu)  +  (Vu)r)  +  Vp  +  A11  •  Vu  =  A f  in  Cl, 

V  •  u  =  0  in  Cl , 

and 

u  =  A(b  +  g)  on  T . 
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where  f  is  a  given  body  force,  b  and  g  are  boundary  velocity  data  with  /r  b  •  n  dT  =  0  and 
J*r  g  •  n  dT  =  0,  and  u  dehotes  the  (constant)  kinematic  viscosity.  We  have  absorbed  the 
constant  density  into  the  pressure  and  the  body  force.  If  the  variables  in  these  equations 
are  nondimensionalized,  then  A  is  simply  the  Reynolds  number  Re.  Given  a  desired  velocity 
field  uo,  we  define  for  any  u  and  g  the  functional 

■7(u,  P,  g)  =  J  Jn  |u  -  uoi4  iO.  +  5  f  (I' v.gp  +  |gp)  dT , 

where  V,  denotes  the  surface  gradient.  We  then  consider  the  following  optimal  control 
problem  associated  with  the  Navier-Stokes  equations: 

min{*7(u,  p,  g ) :  (u,p)  6  2,  g  6  0} 

subject  to  (u,  g)  satisfying  the  Navier-Stokes  equations , 

where  2  and  ©  are  suitable  sets.  For  this  problem,  the  abstract  framework  can  be  used  to 
show  that  optimal  solutions  exist  and  that  the  Lagrange  multiplier  rule  may  be  used,  to 
derive  an  optimality  system  from  which  optimal  states  and  controls  may  be  deduced,  and 
to  derive  optimal  error  estimates  for  conforming  finite  element  approximations  of  solutions 
of  the  optimality  system. 


Feedback  control  of  Karman  vortex  shedding 

The  control  of  the  forces,  e.g.,  lift  and  drag,  exerted  on  a  submerged  obstacle  by 
the  fluid  that  surrounds  it  is  important  in  many  applications.  Here  we  focus  on  the 
problem  of  controlling  the  lift  in  the  plane,  unsteady,  incompressible,  viscous  flow  around 
a  circular  cylinder.  The  control  of  flows  has  a  long  and  rich  history.  Some  of  these 
past  efforts  have  been  devoted  to  conducting  experiments  and  simulations  with  different 
controlling  mechanisms,  without  any  attempt  to  optimally  design  or  to  actively  change 
these  mechanisms.  For  example,  one  may  consult  textbooks  for  detailed  expositions  of 
such  efforts  in  the  area  of  boundary  layer  control.  Of  course,  there  has  been  great  success 
in  the  design  and  application  of  active  controls,  especially  in  the  area  of  aerodynamic 
controls.  These  efforts,  for  the  most  part,  use  simple  models  to  account  for  the  effect  of 
the  flow  on  the  submerged  object. 

Flows  about  a  cylinder  at  even  moderate  values  of  Reynolds  numbers,  e.g.,  Re  >  47, 
exhibit  an  unsymmetric,  periodic  shedding  of  vortices.  As  a  result,  the  lift  force  exerted 
by  the  fluid  on  the  cylinder  is  periodic  in  time.  Here,  our  goal  is  to  show  the  efficacy  of 
a  simple  feedback  control  law  for  the  reduction  of  the  magnitude  of  the  lift.  The  control 
mechanism  used  to  attempt  to  the  reduce  the  size  of  the  lift  oscillations  is  the  injection 
and  suction  of  fluid  through  orifices  on  the  surface  of  the  cylinder.  The  amount  of  fluid 
injected  or  sucked  through  the  orifices  is  determined,  using  a  simple  feedback  law,  from 
the  pressure  “measured"  at  various  stations  on  the  cylinder.  Thus,  in  the  language  of 
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feedback  control,  the  sensor  determines  the  pressure  at  various  stations  on  the  cylinder 
and  the  actuator  injects  of  sucks  fluid  through  orifices  also  on  the  cylinder.  Other  sensing 
and  actuating  mechanisms  could  also  be  used,  e.g.,  the  vorticity  or  stress  components  on 
the  cylinder  for  sensing  and  rotation  or  shape  modification  of  the  cylinder  for  actuating. 

The  type  of  feedback  laws  used  in  our  study  are  based  on  the  observation  that  differ¬ 
ences  in  the  pressure  on  the  top  and  bottom  halves  of  the  cylinder  should  give  an  indication 
of  the  asymmetric  behavior  of  the  lift.  Such  pressure  differences  are  then  used  to  determine 
how  much  fluid  to  inject  or  suck  through  the  orifices.  No  attempt  is  made  to  systematically 
design  an  “optimal”  feedback  law,  i.e.t  one  that  in  some  sense  does  the  best  possible  job  in 
meeting  the  objective.  However,  the  computational  results  we  have  obtained  indicate  that 
the  simple  feedback  law  we  use  here  is  quite  effective  in  reducing  the  size  of  the  oscillations 
in  the  lift.  Details  may  be  found  in  [16]  and  [18]. 

We  have  investigated  the  control  problem  by  numerically  solving  the  2-D  Navier-Stokes 
equations 

~  u  •  Vu  +•  Vp  =  -^-V2u  and  V  •  u  =  0, 
at  Re 

where  u  and  p  are  the  velocity  and  kinematic  pressure  (pressure  divided  by  density), 
respectively.  The  above  equations  are  nondimensionalized  using  the  cylinder  diameter  d, 
the  free  stream  velocity  U,  and  the  kinematic  viscosity  of  the  fluid  u.  The  Reynolds  number 
Re  is  defined  by  Re  =  Ud/v.  Initial  and  boundary  conditions  on  the  cylinder  surface  are 
imposed  on  the  velocity: 

u(x,  0)  =  0  in  the  computational  domain  Cl 
and 

u(x,  t)  =  g  on  the  cylinder  surface  37  and  for  t  >  0 . 

For  the  uncontrolled  flow  about  the  cylinder,  the  cylinder  surface  is  a  solid  wall  so  that 
in  this  case  g  =  0.  When  controls  are  applied,  this  boundary  condition  becomes  inhomo¬ 
geneous  on  the  portions  of  T  covered  by  the  injection  and  suction  orifices,  i.e.,  gpOat 
these  orifices. 

In  order  to  avoid  complications  with  boundary  conditions  at  infinity,  a  specific  finite 
domain  was  defined.  The  origin  of  the  (z,  y )  coordinate  system  is  located  at  the  center  of 
the  cylinder  and  the  cylinder  has  a  unit  diameter.  The  computational  domain  we  use  is 
the  rectangle  — 5  <  z  <  15  and  —  5  <  y  <  5  with  the  cylinder  excluded.  The  geometry 
of  the  domain  is  sketched  in  figure  below.  The  exterior  boundary  of  the  computational 
domain  consists  of  the  four  sides  of  the  rectangle;  again,  see  the  left-hand  figure  below.  At 
the  inflow  boundary  Tj  the  velocity  is  set  to  the  uniform  value  at  infinity,  i.e.,  u  =  1  and 
v  =s  0,  where  u  and  v  denote  the  z  and  y  components  of  the  velocity  vector,  respectively. 
At  the  outflow  T0  a  vanishing  “stress”  boundary  condition  is  imposed.  Specifically,  if 
t  =  — pn  +  (l/Re)du/dn,  then  on  ro  we  set  ti  =  *2  —  0,  where  f*  and  t*  respectively 
denote  the  z  and  y  components  of  t.  The  vector  t  is  not  actually  the  stress  vector; 
however,  it  has  been  found  that  imposing  such  an  outflow  condition  is  often  more  effective 
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than  requiring  the  true  stress  to  vanish.  Moreover,  it  is  a  more  convenient  boundary 
conditions  to  use  with  the  particular  form  of  the  viscous  term  appearing  in  the  Navier- 
Stokes  equations  given  above.  On  the  top  T{  and  bottom  T*  of  the  rectangle  we  impose 
the  mixed  conditions  v  =  0  and  tj,  —  0.  The  boundary  conditions  are  also  indicated  on  the 
sketch  of  the  computational  domain  given  in  left-hand  figure  below. 


For  the  uncontrolled  case,  the  problem  we  wish  to  discretize  is  then  defined  by  the 
above  equations  with  boundary  conditions  on  I\,  Fot  F/,,  and  I\. 

The  spatial  discretization  is  effected  using  the  Taylor-Kood  finite  element  pair  based 
on  a  triangulauion  of  the  computational  fiow  domain.  A  typical  triangulation  is  depicted  in 
right-hand  figure  above;  it  consists  of  1735  triangles  and  3735  velocity  nodes.  The  velocity 
(pressure)  field  is  approximated  using  continuous  piecewise  quadratic  (linear)  polynomials; 
the  same  grid  is  used  for  both  fields. 

The  semi-implicit  single-step  Euler  method  is  employed  for  the  time  discretization. 
Thir  time-stepping  algorithm  is  defined  as  follows:  given  u”“l,  the  pair  (un,pm)  is  ob¬ 
tained  by  solving  the  linear  system 

- - - - -r  (um“l  ■  V)u  and  V  •  um  =  0 , 

where  S  denotes  the  time  step.  Of  course,  (5)  and  (6)  are  supplemented  by  the  appropriate 
boundary  and  initial  conditions.  This  scheme  has  been  previously  analyzed  and  is  found 
to  be  unconditionally  stabe. 

Using  the  computational  scheme  described  above,  we  have  simulated  the  uncontrolled 
problem  at  Re  =  60  and  Re  =  80.  To  break  the  symmetry  in  the  problem  and  trigger  the 
vortex  shedding,  a  perturbation  is  applied  for  a  short  time  interval.  The  transient  solution 
develops  until  a  periodic  state  is  reached  at  approximately  t  =  70  and  t  —  90  for  Re  —  60 
and  Re  =  80,  respectively.  These  phases  of  the  solution  are  perfectly  captured  by  the  lift 
coefficient  Ci  defined  by 
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where  t*  is  the  y  component  of  the  true  stress  vector  and  9  is  the  angle  along  the  surface 
of  the  cylinder  measured  from  the  leading  edge.  The  evolution  in  time  of  Cl  exhibits  an 
oscillation  having  a  period  T  as  7.0  units  in  time  for  Re  =  60  and  T  ss  6.5  for  Re  =  80, 
leading  to  a  Strouhal  number  St  »  0.14  for  Re  =  60  and  St  as  .15  for  Re  =  80;  see 
the  left-hand  figure  below  for  Re  ~  80.  These  values  are  in  agreement  with  experimental 
results  and  other  numerical  simulations. 


We  have  used  a  variety  of  number  and  location  of  blowing  and  suction  orifices  along  the 
back  side  of  the  cylinder.  One  effective  arrangement  is  to  have  two  suction  slots  centered 
at  =23t/32  on  the  back-side  of  the  cylinder  and  a  single  blowing  slot  centered  at  7T.  We 
place  the  sensors  at  symmetric  locations  on  the  front-side  of  the  cylinder  and  we  sense  the 
pressure  at  these  locations.  The  feedback  law  is  chosen  as  follows: 

u  *  min  -  ?(9 2)1, 0}  s(0)  • 

Here,  ct  and  $  are  constants,  with , 3  being  used  to  limit  the  3ize  of  allowable  controls;  g  is 
a  given  velocity  profile  at  the  orifices.  We  have  performed  computational  simulations  for 
numerous  cases  but  report  only  the  case  of  a  =  30  and  0  =  2.0.  The  control  is  turned  on 
at  t  =  100.  The  results  are  shown  in  right-hand  figure  above  ffom  which  one  ran  conclude 
that  the  simple  feedback  law  can  be  quite  effective  in  reducing  the  size  of  the  oscillations 
in  the  lift. 
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